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Large empty regions of space, called voids, represent the 
majority of the volume of the prese nt Universe. Tliey were 
first d is covered in ob s ervati ons by Gregorv fc ThompsonI 
(|l978h . IJoeveer et all (Il978h an d ITuUv fc FisheJ 
followed later bv iKirshner et al] (Il98ll) and more largel y 
in the Cf A redshift catalogue (|de Lapparent et al.l Il986l ) . 
This discovery was followed by a large amount of theo- 
retical w ork. The first gravita t ional instability model was 
given bv iHoffman fc ShahamI j 19821 ). quickly followed by 
iHoffman et all (Il983h for an infinite size regular mesh of 
void and bv iHausman et al.l (|l983l ) for the impact of cos- 
mology on their evolution. Other work studied the general 
self-similar evolution of v oids in Einstein-de-Sitter universes 
(|Bertschingedll983l . [l985l ). 

However, as the voids are intrinsically large and the 
surveys at that time were small, we only detected a small 
number of them. This has hindered their use as a cosmo- 
logical probe for a long time [except some constraints on 
thei r maximal size compatib le with CMB observations by 
e.g. iBlumenthal et all (| 19921 ')]. This situation has changed 
with the advent of deep and wide galaxy surveys s uch as 
the Sl oan Digital sky sur vey (SP SS. lYork et al.ll2000l ). 2dF- 
GRS (ICoUess et al.ll2001^ 2MRS (|Huchral2000l ) and now the 



ABSTRACT 

Wc propose a new definition of cosmic voids based on methods of Lagrangian orbit 
reconstruction as well as an algorithm to find them in actual data called DIVA. Our 
technique is intended to yield results which can be modeled sufficiently accurately to 
create a new probe of precision cosmology. We then develop an analytical model of 
the ellipticity of voids found by our method based on Zel'dovich approximation. We 
measure in A^-body simulation that this model is precise at the ~0.1% level for the 
mean ellipticity of voids of size greater than ^4ft,^^Mpc. We estimate that at this 
scale, we are able to predict the ellipticity with an accuracy of ^ 0.02. Finally, 
we compare the distribution of void shapes in A''-body simulation for two different 
equations of state w of the dark energy. Wc conclude that our method is far more 
accurate than Eulcrian methods and is therefore promising as a precision probe of 
dark energy phenomenology. 

During the last decade, several algorithms to find voids 
have been built. They are separated in three broad classes. 
In the first class, the void finders try to find regions empty 
of galaxies ([Ka uffmann fc F airalll Il99l l: 'El-Ad ct al."l997|; 

NclsonI 



Hovle fc Vogelev 2002 ; Patiri et all 12006 : .Foster fc 
20091 ). The second class of void finders try to identify voids 



as geometrical str uctures in the dark matter distribution 
traced by galaxies JPlionis fc Basilakoj 20021: Colberg et al 



20051 : IShandarin et all l2006l : iPlaten et all l2007l : I Nevrincld 



20081 ). The third class identifies structures dynamically by 



checking gravitati onally unsta b le po i nts in the distribution 
of d ark matter (|Hahn et al.l l2007l : iForero-Romero et al.l 
l2009l ). At the same time, A'^-body simulations focused on 
the stu dies of voids in a cosmologi c al context were flour- 
ishing iMartel fc WassermanI [1990; iRegos fc Gelled [l99ll: 



van de Weygaert fc van KampenI Il993l: J Gottlober et al.l 
20031: 'Benson et al.| l2003l : IColberg et al.1 |2005| ). Recently, 



6dFGS (| Jones et al.ll2009l ). Still we miss a clear and simple 



definition of voids that would allow us to use them as a pre- 
cision cosmological probe. In this paper, we investigate, an- 
alytically and numerically using A^-body simulations, a new 
algorithm for finding voids in Large Scale structure surveys 
and an analytical model that accurately predicts the proper- 
ties of voids found by this method as a function of cosmology. 



Colberg ct al. (200g) made a comparison which shows that, 
even if these currently available void finder techniques find 
approximately the same voids, the details of the shapes and 
sizes found by each of the void finders may be significatively 
different. This problem is further enhanced by the existence 
of ad-hoc parameters in most of the existing void finders, 
which changes the exact definition of voids and does not 
allow reliable cosmological predictions. One aspect of voids 
that is also often left aside is the hierar chical structure s 
of voids. So far, apart from ZOBOV (|Nevrinckl |2008| ) 
and the related W atershed Void Finder (WVF) method 
iPlaten et al.1 |2007| ). which are parameter free, no void 
finder tries to identify correctl y the hierarchy of voids- 
in-voids and clouds-in-voids (|Sheth fc van de Weygaert 
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|2004 ). Another problem of these void finders comes from 
their Eulerian nature; they try to find structures that are 
not necessarily in the same dynamical regime (linear or 
non-linear), which complicates the building of an analytical 
model. 

We propose studying a new void finder that belongs 
to the third class of these void finders. It is based on the 
success of both the Monge-Ampere-Kantoro vitch (MAK) 
reconstruction of the orbi t s of galaxies jBrenier et al.l 
I2OO3I : iMohavaee et al.l I2OO6I: iLavaux et all [ioosl ) and the 
Zel'dovich approximation ( Zel'dovichl Il970l ). This method 
is based on finding a way to compute the Lagrangian coor- 
dinates of the objects at their present position. The study of 
voids in Lagrangian coordinates is not new. The evolution 
of voids in the adh esion approximation has been studied by 
ISahni et all \l994 ) to understand the formation and evolu- 
tion of voids and t heir inner substructure in a cosmolog- 
ical context. Later, ISahni fc ShandarinI (|l996l ) emphasized 
the precision of the Zel'dovich approximation for studying 
void dynamics compared to higher order perturbation the- 
ory, either Lagrangian or Eulerian. However, no void finder 
method has yet tried to take advantage of the Zel'dovich ap- 
proximation for detecting and studying voids in real data. 
The voids detected with this method are going to be intrin- 
sically different than the one found using standard Eulerian 
void finder. This hardens the possibility of making a void- 
by-void comparison of the different methods. 

The use of Lagrangian coordinates gives one imme- 
diate advantages compared to standard void finding: the 
Lagrangian displacement field is still largely in the linear 
regime even at 2 = and especially for voids. This allows us 
for the first time to make nearly exact analytical computa- 
tion on the dynamical and geometrical properties of voids in 
Large scale structures. The MAK reconstruction is thus par- 
ticularly adapted to study the dynamics of voids. However, 
there is an apparent cost to pay: we lose the intuitive way of 
defining voids as "holes" in a distribution of galaxy, that is 
the place where matter is not anymore. On the other hand, 
we gain the physical understanding that voids correspond to 
regions from which matter is escaping. 

The dynamics of voids may provide a wealth of infor- 
mation on dark energy without the need for any new survey. 
The first obvious probe of dark energy properties comes from 
the study of the linear growth factor. Its evolution with red- 
shift depends, among the other cosmological parameters, on 
the equation of state w of the Dark Energy. In this work, 
we assume that w is independent of the redshift. We note 
that in galax;y surveys, our method is going to be sensitive 
to bias but not more than the direct approach to void find- 
ing. Indeed, void finders of the first class are sensitive to the 
selection function of galaxies. Generally this is done by limit- 
ing the survey to galaxies with an apparent magnitude below 
some designated threshold. Changing this selection function 
of the galaxies acts on the boundaries of the detected voids, 
which thus changes the geometry of these voids. From the 
point of view of void finders, this will also act as a "bias". 
The method that we propose has a more conventional de- 
pendence on the bias by using the dark matter distribution 
inferred from the galaxy distribution. The advantage is that 
this bias could be calibrated. One exact calibration consists 
in comparing peculi ar velocities reconst ructed using MAK to 
observed velocities (|Lavaux et al.|[2008l ). Additionally, there 



are a number of other complementary ways of determinin; 



bias from galaxy redshift surveys (e. g. Benoist et al. 199' 



Norbcrg ct al l I2OOII: iTeg mark ct al. ' '2004; [ Erdogdu et al.l 



2005,,; .Tcgmark et al.|[2006. ; Pcrcival ct al...2007l '). 

This paper is a first of a series studying the proper- 
ties of voids found by our void finder. It is organised as 
follows. First, we recall the theory of the Monge-Ampere- 
Kantorovitch reconstruction in Section [2] Then, we explain 
how we can use reconstructed orbits as an alternative way to 
detect and characterise voids. This corresponds to the core 
of DIVA, our void finder through Dynamical Void Analysis, 
and is explained in Section O In Section 21 we model an- 
alytically the voids found by DIVA. In Section [5] we test 
our void finder on Ai'-body simulations. We also check our 
analytical model against the results of the simulations for 
two cosmologies. In Section [6] we compare DIVA to earlier 
existing void finders. In Section [71 we conclude. 



2 THE MONGE-AMPERE-KANTOROVITCH 
RECONSTRUCTION 

The Monge-Ampere-Kantorovitch reconstruction (MAK) is 
a method capable of tracing the trajectories of galax- 
ies back in time using an approximation of the complete 
non-linear dynamics . It is a Lagrangian method, as PIZA 
Croft fc Gaztanagal 1 19971 ') or the Least-Action method 
PeeblesI ligSgl T The MAK reconstruction is discussed in 



grea ,t detail inlBrenier et all (|2003l '). iMohavaee et al] (|2006l ') 
and lLavaux et a" ( 2008l l. Itis based on the hypothesis that, 
expressed in comoving Lagrangian coordinates, the displace- 
ment field of the dark matter particles is convex and poten- 
tial. Since then, this hypothesis has been justified by the 
success of the method on Ai'-body simulations. With the lo- 
cal mass conservation, this hypothesis leads to the equation 
of Monge- Ampere: 



deti 



9^$ 
' dqidqj 



P(x(q)) ^ 
Po 



(1) 



with q the comoving Lagrangian coordinates, x(q) the 
change of variable between Eulerian (x) and Lagrangian co- 
ordinates (q), p(x) the Eulerian dark matter density and po 
the initi al comoving density o f the Universe, assumed homo- 
geneous. |Brenier|eFal] (j2003h showed that solving this equa- 
tion is equivalent to solving a Monge-Kantorovitch equation, 
where we seek to minimise 



/[q(x) 



dV(x) (x-q(x))^ 



(2) 



according to the change of variable q(x). Discretising this 
integral, we obtain 

\ 2 



<T(i); 



(3) 



with (T a permutation of the particles, qj distributed ho- 
mogeneously, Xi distributed according to the distribution of 
dark matter. Doing so, we obtain a discretised version of the 
mapping q — >■ x(q) on a grid. 

To solve for the problem of minimising Eq. ([Sj with 
respect to a, we wrote a high-performance algorithm that 
has been parallelised using MPI. T his algorithm is b ased on 
the Auction algorithm developed bv lBertsekad (| 19791 ). It has 
an overall time complexity for solving cosmological problems 
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empirically between 0{n^) and 0{n'^) (at worst) with n the 
number density of mesh elements {qj}u 



3 THE VOID FINDER BY ORBIT 
RECONSTRUCTION 

In this section, we describe our void finder DIVA (for Dy- 
namical Void Analysis). First, we define in Section [3.11 what 
we call a void in this work. Second, in Section [3.21 we make 
use of the displacement field in the immediate neighbour- 
hood of a void to define the ellipticity arising from tidal field 
effects, which we also call tidal ellipticity. In Section [3.3l we 
define the Eulerian ellipticity of our voids. In Section [3.4l we 
discuss the impact of smoothing in Lagrangian coordinates 
to compute void properties. 

In later sections, we use pure dark matter A''-body simu- 
lations to check the adequacy of the voids found using MAK 
reconstructed displacement field and the one detected in the 
simulated displacement field. The results given by the an- 
alytical models are then compared to the one given by the 
simulated field for two equations of state of the Dark Energy. 

3.1 Definition of a Void 

So far, voids have only been described using a purely geo- 
metrical Eulerian approach. Typically, as mentioned in the 
introduction, a void is an empty region delimited by either 
sphere or ellipsoid fitting or by using isodensity contours. 
We propose here to use a Lagrangian approach and use the 
mapping between Lagrangian, q, and Eulerian coordinates, 
X as a better probe for voids. In the rest of this article, 
we will consider these two coordinates to be linked by the 
displacement field 

x(q) = q-f *(q). (4) 
We now define the source of the displacement field by 

As the displacement field is taken to be potential it is strictly 
sufficient to only look at to study ^. 

We now define the position of a candidate void centre by 
looking at maxima of 5* in Lagrangian coordinateslf] This 
will effectively catch the source of displacement and from 
where the void is expanding. The other, practical, advan- 
tage is that Sis is quite close to the opposite of the linearly 
extrapolated initi al density perturbation s of the considered 
patch of universe (jMohavaee et al.ll2006l ). So we can use the 
usual power spectrum to study most of the statistics of this 
field. So the main approximation we use in the rest of this 
study is that the primordial density field power spectrum 
is a good proxy for the power spectrum of the seed of dis- 
placement and that this displacement is a Gaussian random 
field. 

^ An implementation of the algorithm is cur- 
rently available on the first author's website 
I htt p : / / www . iap . fr /users /lavaux /code . php [ 

^ This maxima corresponds in terms of the primordial density 
field to what is sometimes called a protovoid ijBlumenthal et al] 
I1992I; iPiran et al.]|l993l: iGoldwirth et al.lll995h . 



From we define the matrix Ti^rn. of the shear of the 
displacement, which is linked to the Jacobian matrix: 

dxi 9*; 

Ji.miq) = t; = 5i,m + t; — (q) = 5l,m + Ti.m , (6) 

oqm. Oq,n 

^(q) = \Jl,m\, (7) 
with 

r,„(q)^£i. (8) 

J is the Jacobian of the coordinate transformation q — > x. 
Geometrically, J specifies how an infinitely small patch of 
the Universe expanded, in comoving coordinates, from high 
redshift to z; = 0. We put Ai(q) the three eigenvalues of 
21,m(q) and sort them such that Ai(q) > A2(q) > A3(q). 
Among the candidate voids, we select only voids that have 
strictly expanded, which equivalently means that J > I. 
We may now define three classes of voids that are inspired 
from the usual classes of observable large scale structures 
for galaxies: 

- true voids for which Ai > 0, A2 > 0, A3 > 0. These 
should be the most evident and easily detectable voids as 
they consist in regions which are expanding in the three 
directions of space. 

- pancake voids for which Ai > 0,A2 > 0,A3 < 0. The 
pancake voids are closing along one direction of space but 
expanding along the two other directions. With a geomet- 
rical analysis, this case cannot be distinguished from the 
true void case. However the dynamical analysis is capable of 
that, and this will cause a crucial difference to the analysis 
as we will see later. In practice they represent a substantial 
fraction of the voids. 

- filament voids for which Ai > 0, A2 < 0, A3 < 0. 

We refer to Section 14.31 for the quantitative relative num- 
ber of voids in each of class. As we will see, the distinction 
between those cases is important to quantify the shape and 
properties of voids that we observe at the present time. We 
discuss our definition of voids, and compare it to other void 
finders, in Section [6] 

Here, we have not yet touched the issue of defining the 
boundary of a void. The exact study of the properties of 
void volumes will be studied in a forthcoming paper (Lavaux 
& Wandelt 2009, in p reparation). We us e a variant of the 
Watershed transform (jPlaten et al.|[2007l ) to define the La- 
grangian volume of a void. In Lagrangian coordinates, the 
voids occupy half of the volume. So, instead of enforcing 
strictly that we should have complete segmentation of the 
volume in terms of void patches, we impose that voids must 
correspond only to the places that are sources of displace- 
ment and not to sink like clusters. Gontrarily to a pure Wa- 
tershed algorithm, we thus enforce that > everywhere 
within the void boundary. 

3.2 Ellipticity of a void 

After having defined the position and the dynamical prop- 
erties of the void, we may defi ne an interesting property of 
those structures: the ellipticity. Ilckd (|f 9841 ) first emphasized 
that isolated voids should evolve to a spherical geometry. 
But, in the real case, voids are subject to tidal effects. As- 
suming the present matter distribution evolved from a to- 
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Initial conditions Present time 

Figure 1. Picture of a void and the formation of intrinsic ellipticity - We represent in this figure tlie central idea of tlie definition 
of a void and its ellipticity. We take voids as maxima in S^s,. They correspond to first order to minima of the primordial density field 
represented here by painted surface. These minima undergo an overall expansion from initial conditions to present time. The shape of 
the Void is defined locally at the minimum. The ellipticity is equal to the square root of the ratio of the axes of the ellipsoid which locally 
fits the surface. 



tally isotro pic and homogene ous distribution, IPark fc Led 
lj2007h and iLee fc ParkI (j2009l ) have shown that the distri- 
bution of the ellipticity which is produced by tidal effects is 
a promising probe for cosmology. More generally, previous 
work have shown that a lot of potentially observable statis- 
tical prop erties of voids are directly to the primordial tidal 
field (e.g. iLee fc ParkI l2006l : IPlaten et al1l2008l : IPark fc Leel 
|2009b.a). However , questions may be raised by the direct 
use of the formula iDoroshkevichI (|l970l ') , as they applied it 
in Millennium Simulation. Using the orbit reconstruction 
procedure, our approach should be able to treat the prob- 
lem right from the beginning, even in redshift space (see 
iLavaux et al.|[200^ for a long discussion), though some care 
should be taken for the distortions along the line-of-sight. 

The other advantage of Lagrangian orbit reconstruction 
is that it offers for free a way of evaluating the ellipticity 
locally, potentially at any space point. From the mass con- 
servation equation and the definition of eigenvalues of Ji^m 
we may write the local Eulerian mass density as: 



Peiq.) = 



P ^ pVl 

1(1 + Ai(q))(l + A2(q))(l + A3(q))j VE(q) 



, (9) 



with Vl the Lagrangian volume of the cell at q and VE(q) 
the Eulerian volume of this same cell, p the homogeneous 
Lagrangian mass density. This equation is valid at all times. 
Now we may also explicitly write the change of volume of 
some infinitely small patch of some universe 



1/E(q) = Vi 1(1 + Ai(q))(l + A2(q))(l + A3(q))! 



(10) 



Provided the eigenvalues Xi are greater than -1, which is 
always the case for voids, we may drop the absolute value 
functionlfl Now, with the analogy of the volume of an ellip- 



An eigenvalue less than -1 would mean that the void would 
have suffered shell crossing at the position of its centre, which is 
dynamically impossible as we are at the farthest distance possible 
of any high density structure. 



soidQ we may write the ratio v between the minor axis and 
the major axis 



1 + A3(q) 



(11) 



1 + Ai(q) ' 

and the ratio fj, between the second major axis and the major 



M(q) 



1 + A2(q) 



1 + Ai(q) ■ 
This allows us to define the ellipticity 



e(q) = 1 - i.(q) = 1 - 



1 + A3(q) 
1 + Ai(q) 



(12) 



(13) 



We will define the ellipticity of a void as the value taken by 
e at the Lagrangian position of the void. 

A picture of the concept of voids and ellipticities in this 
work is given by Fig. [T] The painted paraboloid represents 
a small piece of a larger 2D density field whose value is 
encoded in the height and the colour. According to our def- 
inition, the void is at the centre of the paraboloid. At this 
centre, the surface of the volume element is mostly circular. 
The tidal forces are locally transforming the shape of this 
surface, which produces the new elliptic shape on the right 
side of the figure. The surface has been here extended along 
one direction and slightly compressed along the other. 

Though we are not strictly limited to study elliptic- 
ity only at the position of the void it may be promising 
in terms of robustness to non-linear effects. Indeed, due to 
the absence of shell-crossings inside voi ds, the MAK recon- 
struction should give the exact solution terenier et al.ll200^ : 
iMohavaee et al.l l2006l : iLavaux et al] |2008| ) to the orbit re- 
construction problem. It means that the ellipticity that we 



* We used here the convention of iPark fc Led ll2007l'l who take 
the square root of the ratio to define fi and u. 
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will compute will be exact, to the extent that we have taken 
care of the other potential systematics due to observational 
effects (jLavaux et alj [20081 ) . As any other method relying 
on dark matter distribution, we will be sensitive to the fact 
that the large-scale galaxy distribution is potentially biased. 
However, if the bias does not depends wildly on redshift, we 
should be able to compute statistics on ellipticities and de- 
rive the evolution of the growth factor of Large Scale Struc- 
tures. 

We note that, using MAK reconstruction, we have ac- 
cess to the joint distribution of the three eigenvalues. Our 
computation of the ellipticity consists in a projection of the 
whole joint 3d joint distribution on a Id variable. For cos- 
mological analysis, it is not entirely clear which estimator is 
the more robust. On one hand, our intrinsic variables are the 
eigenvalues and we could include them in the analysis just 
as well as the ellipticities On the other hand, using ellipticity 
may be helpful to average on a lot of different voids. It may 
be a more robust estimator with respect to badly modelled 
tails of the distribution of eigenvalues. In this work, we focus 
on the use of the ellipticity, as defined in Ea. (|13p . 

3.3 Eulerian ellipticity 

We define the volume elli pticity Syoi using the e igenvalues of 



the inertial mass tensor (jShandarin et al.ll2006l l: 



(14) 



Mxx = ^ mi{y1 -f zf), Mxy = — ^ irnxiyi, 

1=1 1=1 
where rrii and Xi, yi, Zi are the mass and the coordinates 
of the i-th particle of the void with respect to its centre 
of mass. The other matrix elements are obtained by cyclic 
permutation of x,y and z symbols. We put Ij the j-th eigen- 
values of the tensor M, with 7i ^ I2 ^ Is- We may now 
define the volume ellipticity as in 



h + h- h 
I2 + I1- h 



1/4 



(15) 



Even though our work is focused on the tidal ellipticity (Sec- 
tion [5]2J, there is some interest to compare how the Eulerian 
volume ellipticity compares to the local tidal ellipticity, as 
most of the existing void finders use EvoI as a probe for the 
dynamics. 

To have a fair comparison with DIVA results, we are 
computing the inertial mass tensor from the displacement 
field ^(q) smoothed on the scale as for the rest of the anal- 
ysis. The void domain is defined as specified in Sec. 13. II The 
inertial mass tensor is thus: 



I d'^q {{qy + *,(q))2 + (q, + vl/.(q))^) , (16) 



d'q(g. + *4q))(g, + *,(q)). (17) 



with the other elements obtained by cyclic permutations. 
The volume ellipticity Svoi is compared to the tidal elliptic- 
ity emvA in Section 15.3.31 Except in that section, we only 
consider sdiva in this paper. 

3.4 Smoothing scales 

There is an apparent price to pay to go to Lagrangian co- 
ordinates. One has to set a smoothing scale in Lagrangian 



coordinates and study the dynamics at corresponding mass 
scale and let go of the evident notion of whether we see a 
hole in the distribution of galaxies or not. It actually could 
be an advantage. Smoothing at different Lagrangian scales 
allows to probe the structures at different dynamical epoch 
of the void formation. Each Lagrangian smoothing scale cor- 
responds to a different collapse time: the smallest scales be- 
ing the fastest to evolve. DIVA in this respect allows us to 
study the dynamical properties of a the voids which have 
the same collapse time. 



p eak p a 

of structure formation (IBond fc Mverj Il996l '). which is a 
simplified but quite accurate model of the dynamic of 
peaks in the density field. This model is even more pre- 
cise for the void patches, which is the n ame of the equiv- 
alent model for studying voids ( see e.g. Sahni ct al.l Il994l : 
ISheth fc van de Wevgae"3 |2004| : iNovikov et al. 20(3). Of 
course, the number of voids depends on the filtering scale 
(see Section 14.31 and Section I5.3.2|l . If we smooth on large 
scales we should erase the smaller voids and leave only the 
voids whose size is large enough. 

Smoothing also affects the ellipticity distribution. As 
we smooth to larger and larger scales the density distri- 
bution probed by the filter should become more and more 
isotropic. This leads voids to become more spherical and 
thus the ellipticity distribution should be pushed towards a 
perfect sphere. In this paper, we consider a few scales sep- 
arately and try to understand what were the properties of 
the minima at each of these scales (see Section [5. ip . 



4 ANALYTICAL MODELS FOR VOIDS 

In this section, we describe an analytical model of the dis- 
placeme nt field. This mod el is based on Zel'dovich approx- 
imation (IZel'dovichl 1 1970l ) . In a first step (Section 14. ip . we 
recall the statistics of the shear of the displacement field. 
Then, in Section 14.21 we express the ellipticity defined by 
Eq. (|13p in terms of this statistic. Finally, we explicitly write 
the required statistical quantity in the model of Gaussian 
random fields and give some expected general properties of 
the voids in this model in Section [4.31 



4.1 Displacement field statistics 

iPark fc Led (|2007h described an analytical model of void 
ellipticities based on the Zel'dovich approximation. This 
model should be particularly suitable on making predictions 
of the result given by D IVA, knowing the previous successes 
of MA K in this domain (iMohavaee et al.| | 2006l : lLavaux et al.l 
12008}). The model that IPark fc Led (|2007l 'iliave proposed is 
based on the unconditional joint distributio n of the ei genval- 
ues of the tidal field matrix J;_m (Doroshkcvichll 1970l ) . given 
the variance of the density field (Appendix E)) : 



P(Ai,A2,A3la) = 



3375 



exp 



3 (2Kf - 5A"2 



2(72 



!(Ai-A2)(Ai-A3)(A2-A3)|, 



(18) 
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with 
K2 



^l^'S + -^l'^2 + A2A3 . 



(19) 
(20) 



as defined in Appendix [XI 

This expression however neglects the fact that voids cor- 
respond to maxima of the source of displacement [f] As the 
curvature of 5* = Ai + A2 + A3 is correlated with Ji,m, we 
need to enforce that we are actually observing the eigen- 
values in regions where the curvature of 5« is negative. A 
better expression would be derived if we could constrain 
that the Hessian H (the matrix of the second derivatives) 
of 5** is negative, which is the case in the vicinity of max- 
ima of Sijj, the source of the displacement field. We derive 
in Appendix |B] a general formalism that allows us to com- 
pute numerically the probability P(Ai, A2, A3|(Jt, r, _ff < 0) 
to observe the eigenvalues {Ai, A2, A3} given that we look in 
these regio ns. This form a lism i s a natural extension of the 
formula of iDoroshkevichI (|l97Cll ) (for which a simple deriva- 
tion is given in Appendix [X)) . 

Of course, "true voids" have the additional constraint 
that Ai > for all i = 1,2,3. As we assumed in previous 
sections that eigenvalues are ordered according to Ai > A2 > 
A3, the constraint A3 > is sufficient to study this case. 

4.2 Ellipticity statistics 

Whether we use the conditional probability 
P{Xi, X2, X3\aT,r, H < 0) or the unconditional one 
-P(Ai, A2, Aalor), both under notation ^(Ai, A2, Aalor, r), 
we may now express the probability to observe S, v, jj, 
[defined in Equations ((SJ, (fTTj) and (fT^ j in terms of V: 

P{fi, u, 5\r, ar) = P(Ai, A2, Asjr, ar) x 4(<? - 3) V'^ ^ (gl) 



with 

Ai = 

A2 = 
A3 = 



1 - 2/1^ + 5fi^ + ly^ 
-2 + 5 + + 1.2 



The ellipticity distribution of voids is thus 



(22) 
(23) 
(24) 



-1 r + oo pi 

P{e\oT,T) = — / P{^,,u,5\aT,r)A^iAS, (25) 



with 



f dfidud5P{fi,v,S\aT,r) . (26) 

(5— — 00 J ly— J ^—v 

The alternative distribution for "true voids" is given by en- 
forcing that Ai > and may be obtained by introducing the 
Heaviside function 0(Ai) in Eq. (|2ip and renormalising. 

We note that the ellipticity that we are considering 
here is of dynamical nature (as emphasized bv lPark fc Led 



In terms of primordial density fluctuations, voids correspond to 
minima of the density field. As MAK is providing a good approx- 
imation of this field, we may safely jump from one concept to the 
other. 



This comes in contrast with the first studies of void 
ellipticities due to reds hift distortions bv lRvdenl(|l995h and 
iRvden fc Melot^ (|l996t ). We do not discuss this earlier defi- 
nition of ellipticity but only the later one. 

4.3 Application to cosmology 

Shapes of voids ~ Now we may compute the ellipticity distri- 
bution of voids, P{e\S\s>) for a given cosmology. The variance 
of the density field ax, assuming the power spectrum of the 
density fluctuations P{k), is given by 



2 

CTx = 



with 



1 r + oc 

P{k)Wl^ik)k^dk, 



WRjik) =exp ( 



(27) 



(28) 



the Fourier transform of the filter function used to smooth 
the density field on the scale Rf. In practice, we smooth 
the displacement field in Lagrangian coordinates, to reduce 
noise and non-linear effects appearing at small scales in the 
MAK reconstruction (< 5/i~^Mpc). Thus, we will compute 
the ellipticity distribution of voids, given that we smoothed 
on the scale Rf in Lagrangian coordinates, and that the 
local source of displacement of the void is ^^(q). 

With the model P(Ai, A2, AsIctt, r, < 0) of Ap- 
pendix[Bj we may also estimate the number of voids in each 
class we defined in Section [3]T] As an illustration, we picked 
a usual ACDM cosmology, with Q.^ = 0.30, fib = 0.04, 
(jg = 0.77, h = 0.65 and estimated the fraction of voids in 
each class. The results are, when we smooth at 4/i~^Mpc: 

- true voids: We estimate that these voids represent 
~40% of the primordial voids, which correspond to under- 
densities in the primordial density fluctuations. 

- pancake voids: Doing the same estimation as for "true 
voids" , we get that ~50% of the primordial voids should be 
in that case. 

- filament voids: They correspond to ~10% of the pri- 
mordial voids. 

We show in Figure [2] the analytical distributions 
of ellipticity for the two cases when one constrains (or 
not) the curvature of to be negative. The solid curve 
corresponds to the ellipticity distribution obtained using 
P{Xi, A2, XalH < 0). The dashed curve is obtained with the 
unconstrained distribution. In the left panel, we took all 
voids with Ai > —1. In the right panel, we restricted our- 
selves to "true voids". The difference is striking in the left 
panel between the two models, whereas in the right panel 
they essentially give the same prediction. This can be un- 
derstood by looking at the value of the correlation coeffi- 
cient between the curvature of and T;,„i (also defined in 
Eq. dmoj) 

(29) 



V 5256 



with 



1 r+oc 



(30) 



This coefficient is equal to ~-0.67 for the aforementioned 
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Figure 2. Comparison of MAK reconstructed and analytical ellipticity distribution - We represent here the distribution of eUipticity of 
the voids, marginalised over all possible 5\p. We used black squares for the ellipticity distribution obtained using the MAK reco nstructed 
displacement field applied on the simulation. The dashed blue curve is computed using the unconditional IPoroshkevicbl |l97^) formula. 
The red curve is our new formula obtained by conditioning that voids are regions where the density field has a positive curvature. The 
left panel gives the result for all voids (true, pancake and filament). The right panel gives the result for only true voids. The blue dashed 
curve and red solid curve overlap. All fields were smoothed with a Gaussian kernel of radius 4/i~^Mpc. 



cosmology. This indicates that the two curvatures are quite 
strongly linked. Thus, enforcing that T;,™ is positive causes 
the curvature of to be preferentially negative. So, the two 
distributions of the right panel of Fig. [5] should be similar. 
On the other hand, for the left panel no such implication 
exists, which leads to the largely evident discrepancies of 
ellipticities. We note the distributions of the right panel are 
only measurable using our void finder as the other ones can- 
not distinguish truly expanding voids and pancake voids just 
by looking at their shape. 

Number of voids - Now that we know the shape the 
voids should have in the context of linear theory, we would 
like to know how many of them should be present in the Uni- 
verse. With our definiti on of voids, we may co nveniently use 
the results obtained by iBardeen et al.l (|l986l ) for Gaussian 
random fields. In particular in the equation (4.11b), they 
show that the average number density of maxima is equal 
to 



29-6^ 
53/22(27r)2i?3 



: 0.016i?: 



with 



(31) 



(32) 



in our notation. We note that this is a mean number, so 
we expect some fluctuation according to the mean which 
are difficult to compute analytically. Tests on Gaussian ran- 
dom field seems to point out that the number of voids are 
distributed like a Poisson distribution. We also expect this 
number to slightly overestimate the actual density of voids 
that we will find in simulations (Section 15. 2|) . 



In the next section, we are now going to confront the 
analytical model with the results given by DIVA applied on 
A''-body simulations. 



5 TEST ON iV-BODY SIMULATIONS 

In this section, we are going to identify voids in the A^-body 
samples described in Section lSTT] We then give a sketch (Sec- 
tion [52]) of the procedure we followed to perform the MAK 
reconstruction, whi ch corresponds to the one described in 
iLavaux et al.l (|2008l ). In Section [5.31 we focus on the results 
obtained at z = 0. First, we give an illustration of a void of 
each class in Section [537l] We then compare the results ob- 
tained using the displacement field given by the simulation 
and the one reconstructed by MAK (Section I5.3.2p . There, 
we also detail the number of voids detected and their ellip- 
ticities for both fields. We compare quantitatively the dis- 
tribution of Eulerian volume ellipticity to the Lagrangian 
tidal ellipticity in Section [5.3.31 Finally, we check the valid- 
ity of the analytical model on the simulated displacement 
field (Section 15. 3. 4p . In Section [531 we look at the evolution 
of the mean ellipticity for a simulation with lo = — 1 and in 
the analytical model. At last, in Section 15.51 we assess the 
possibility of measuring the evolution of the mean ellipticity 
in two simulations where w is either —1 or —0.5. 



5.1 The A'^-body simulations 

To test our void finder, we use three large volume TV-body 
simulations but with medium resolution of A'^ = 512"^, L — 
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SOO/i^^Mpc, fim = 0.30, = 0.70, = 65 km s"^ Mpc-\ 
ris = 1, (78 = 0.77, f2b = 0.04. The A''-body simulations con- 
tains only dark matter and we include the eflFect of baryons 
only through power spectrum features in the initial condi- 
tions. This essentially reduces power on scales smaller than 
the sound horizon and introduces Baryonic Acoustic Oscilla- 
tions. The two first A^-body simulation (ASIM and ASIM2) 
corresponds to a standard ACDM cosmology for which the 
equation of state is given by w — —1. ASIM and ASIM2 
have exactly the same cosmology but have different initial 
conditions. They will allow us to assess the impact of look- 
ing at two different realisations of the initial conditions. The 
third, wSIM, is assuming an equation of state w = —0.5 for 
the Dark Energ y. To build the ini tial conditions, we modi- 
fied the GRAFIC ( Bertschingeij[200ll ) packa ge to use the powe r 
spectrum generated by the CAMB package (|Lewis et al.ll2000l ). 
We reach a sub-Mpc resolution scale which is sufficient for 
proper description of most voids (f-20 /i~^Mpc). The in- 
termediately large volume allows accounting for large-scale 
tidal field effects and cosmic variance effects. We used the 
paral lelised version of the RAMSES A-body code (|Tevssieil 
|2002| ) and run it both on the Cobalt NCSA supercluster 
and th e Teragrid NCSA s upercluster through Teragrid fa- 
cilities (|Catlett et al.ll2008l ). We modified RAMSES to simulate 
cosmologies with a dark energy equation of state difi'erent 
than w = —1. 

To account for the impact of clustering, we analyse the 
displacement field for which the mass of dark matter halos 
has been put to the centre of mass of these halos. To be able 
to do that, we construct a halo catalogue using a Friend-of- 
Friend algorithm with a traditional linki ng value of I = 0.2 
iDavis et al.ll 19851 : lEfstathiou et al.ll 19881 ) . We put a thresh- 
old of Ap = 8. This prescription in practice should mostly 
erase the information contained in halos while retaining the 
dynamics outside of them. Though the use of such a low 
number for the particles of halos are questionable for the 
study of the properties of those halos, we are not here in- 
terested in them. We are only interested in checking that 
we keep most of the information useful to study voids and 
their overall dynamics, even though we destroy the small 
scale information. The above prescription has already been 
success fully applied for the study of peculiar velocities with 
MAK (jLavaux et al.ll2008l '). We note that we will only use 
the halo catalogue to do the MAK reconstruction. All the 
tests of the displacement field of the simulation are run on 
the particles of the simulation. 

We note that that the initial conditions of the simu- 
lation present two particularities that must be taken into 
account. The largest mode of the simulation box is fciow = 
1.25 X 10~^ h Mpc~^ thus introducing a sharp low-k cut. 
Additionally GRAFIC applies a Manning filter on the initial 
conditions to avoid aliasing. In practice, GRAFIC multiplies 
the cosmological power spectrum by the following filter: 



,(fc) = 







if kAx ^ n 
elsewhere 



(33) 



with Aa; = 0.976/i^^Mpc the Lagrangian grid step size of 
our simulations. These two features must be introduced in 
the power spectrum of Eq. H27p and (|30p to make correct 
predictions for the simulation. In real observational data, 
no such features should be present. 



5.2 MAK reconstruction and void identification 

Among the different tests that we are going to present in 
the following, we have run only one MAK reconstruction for 
the full halo catalogue. We chose a resolution of 256^ mesh 
eleme nts generated following the algorithm of iLavaux et al.l 
(|2008l ). This algorithm consists in splitting a mass in an in- 
teger number of MAK mesh elements, with the constraints 
that the splitting must be fair and the number of mesh el- 
ements is fixed and equal to 256'^. This method works well 
and has, so far, not been prone to biases. Choosing this num- 
ber of elements gives us a resolution of ~2/i~'^Mpc on the 
Lagrangian coordinates of the displacement field. We can- 
not run it at the full simulation resolution due to the high 
CPU-time cost which hinders doing several reconstruction. 
One reconstruction takes ~13,000 accumulated CPU-hours 
on Teragrid cluster at NCSA. However, as the complexity 
grows as 0{N^''^''), increasing the resolution to 512"^ would 
have consumed ~ 10^ CPU-hours. So we limited ourself on 
running the reconstruction at the 256'' resolution, the ex- 
pected worst case for the performance of MAK. At higher 
redshift, the MAK reconstruction converges faster and gives 
a more and more precise displacement field compared to the 
one given by simulation. These two features are both caused 
by the decrease of small scale non-linearities at earlier times. 
We took the halo catalogue built on ASIM at z — and ran 
a reconstruction on it. The other results presented in this 
paper use the displacement field given directly by the simu- 
lation after having checked that the reconstruction is indeed 
sufficient for our purpose. This is the case as we will not 
look at voids smaller than ~l/i~^Mpc scale in Lagrangian 
coordinates. 

We chose two different smoothing scales on which we 
study the displacement field for voids: 2.5/i~^Mpc and 
4/i~^Mpc. Once the displacement field has been smoothed, 
we compute the eigenvalues and the divergences in Fourier 
space. We then locate the maxima in the divergence of the 
displacement field. At each maxima, we compute the ellip- 
ticity e with the help of Eq. (|13|) . where the Xi are taken as 
the eigenvalues of T;,m(q). The displacement shear tensor, 
is computed numerically from '4'(q) in Fourier space; 



(2^) 



d-'kifc™*,(k)e*". 



(34) 



where ^;(k) is the Fourier transform of the displacement 
field in Lagrangian coordinates. All these quantities were 
evaluated using Fast Fourier Transform on the Lagrangian 
grid. 

In summary, the steps are the following: 

- First, we prepare a catalogue for a MAK reconstruction. 
This involves doing fair equal mass-splitting of the objects. 

- Next, we run the MAK reconstruction. 

- After the reconstruction is finished, we put the com- 
puted displacement field given by MAK on the homogeneous 
Lagrangian grid. 

- We smooth this displacement field and compute the 
tidal field Tij in Lagrangian coordinates in Fourier space 
using Ea. (|34|l . 

- Now we compute (q) on the grid using Eq. ((5)1 , which 
corresponds to summing the three eigenvalues of Tij . 

- We look for local maxima in S*(q) using an iterative 
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steepest descent algorithm. This gives us the list of the voids 
in Lagrangian coordinates. 

- Using these coordinates, we now compute the ellipticity 
using the eigenvalues of Tij at the location of the minima. 

- For computing the void boundary, we execute the mod- 
ified Watershed transform of Section [3. II This identifies the 
Lagrangian domain of the voids. The boundary of this do- 
main is then transported using the displacement field to re- 
cover the Eulerian void volume. 

We now look at the results obtained with MAK, the simu- 
lation and the analytical model at z = in the next section. 

5.3 Results at z — 

5.3.1 Example of found voids 

Before going to the statistical study of the local shape ediva 
properties of void found by DIVA, we look at a visual ex- 
ample of each of the void type. We chose a filtering scale of 
4/i~^Mpc in Lagrangian coordinates. We selected one void of 
each class. These three voids have the following properties: 

- The first void is a true void. The eigenvalues of the 
tidal field r;,™(q) (O at the centre are (1.2,0.84,0.69) 
along the three axis. We thus derive the ellipticity e = 
0.13. The volume of the void, in Lagrangian coordinates, 
is Vl = 75240/i~^Mpc'^, which corresponds to an equiva- 
lent spherical volume given by a sphere of radius -Roquiv = 
(3/(47r) * Vf'"^ = 26/i-^Mpc.' 

- The second void is a pancake void. The eigenvalues of 
the tidal field are (1.1,0.11,-0.60), the ellipticity is 0.563 
and the Lagrangian volume is Vl = 1560/i~^ Mpc'^, with 
iJequiv = 7.2/i-^Mpc. 

- The last void is a filament void. Again, the eigenvalues 
of the tidal field are (0.99,-0.24,-0.61), the ellipticity is 
£ = 0.557 and the Lagrangian volume is Vl = 260/i"^ Mpc^, 
with TJequiv = 4.0/1" ^Mpc. 

Those three voids are represented in three dimensions, along 
with the particles of the simulation in the same region, in 
Fig. [3] We note that the true voids is the largest one. We 
expect to observe that effect as true voids expands in three 
directions and thus are more likely to be bigger than pan- 
cake voids and filament voids. For these three cases, the 
surface delimited by DIVA seems to nicely fit to the struc- 
tures located on the boundaries. In the case of the pancake 
and filament voids, we note that the cavity seem to be split 
into two pieces. This splitting is due to the Watershed trans- 
form criterion which isolated two void cavities separated by 
a structure. 

5.3.2 MAK vs Simulation 

We now concentrate on the properties of the voids at z = 0. 
This is the time where the density is the most non-linear 
and the reconstruction is the most difficult and therefore 
represents a worst case example. We take the displacement 
field given by the simulation as the field of reference to study 
voids. Indeed, this field has been obtained by solving com- 
pletely the equation of motion for each particle. In this sec- 
tion, we will first compare the properties of the voids found 
using this field and the MAK reconstructed field. Then, we 
will focus on how it compares to the analytic model. 



Filter 


Predicted 
average number 


Displacement 
field 


Number of 
candidates 




42908 


Simulated 


31002 






Reconstructed 


28397 


4h-iMpc 


11706 


Simulated 


10643 






Reconstructed 


9412 



Table 1. Unconditioned number of voids in ASIM - we give here 
the unconditioned number of voids found within the volume of 
the simulation, (500 h~^Mpc)^. The predictions arc obtained us- 
ing Equation I I31I I applied on the power spectrum of primordial 
density fluctuations multiplied by the Fourier transform of the in- 
dicated fllter in the flrst column. The last column gives the actual 
number of void candidates that wc found using the displacement 
fleld. 



We represent on Fig. |4] the distribution of ellipticities 
obtained using both the reconstructed and the simulated dis- 
placement field. We also give the number of voids found in 
simulated displacement field, in the reconstructed displace- 
ment field and the expected number of maxima according to 
Eq. (|31|l (Table [T]). To allow for a void-by- void comparison, 
we tried to match the voids found using the two displace- 
ment fields. We considered that two voids are the same if the 
distance between the two voids is less than a Lagrangian grid 
cell (d ^ •\/3Zccii, with Zcoii the length of the side of a cell). 

At 2.5/i~^Mpc smoothing, the agreement of the elliptic- 
ity distribution derived from the simulated displacement and 
the MAK reconstructed displacement is very good, though 
the high ellipticity tail seems a little different in the two 
cases. This is actually due to a selection effect which, unfor- 
tunately, is correlated with the ellipticity. If we look at the 
number of voids detected using the two fields (Table [ij , we 
see that MAK is missing about 10% of the voids detected in 
the simulation. The distribution of ellipticity of those voids 
happens to be skewed towards higher ellipticities. Thus it 
seems that we tend to miss the most elliptical voids. This 
behaviour is expected as these voids tend to be pancake 
voids. So they are closing along one direction, and the more 
elliptical they are the faster they are closing. If they close, 
the particles of these voids shell cross and MAK is not able 
to reconstruct the displacement field. Looking at this same 
distribution, but for a 4h,^^Mpc smoothing, we now barely 
see a difference between the two curves. We indeed checked 
that the ellipticity distribution of the voids that have not 
been identified using the MAK reconstructed displacement 
field is similar to the distribution of the other voids. 

The number of voids detected in the simulation looks 
less than the expected average number of minima (Table [l|. 
This is also due to the destruction of minima by the collapse 
of large scale structures. When we look at the beginning of 
the simulation we find 11,485 minima (for Rf = 4/1" ^Mpc), 
and this number steadily decreases to the value we put in the 
table as the simulation evolves. We estimated the mean and 
the variance in the number of minima using 20 realisations 
of a Gaussian random field with the same cosmology as the 
simulation. We found that the mean should be 11,762 and 
the standard deviation is 58, which is in agreement with the 
result given by the analytic computation. 

Using the match between voids coming from the two 
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fields, we can build a statistical error model in the form of 
the joint probability distribution P(£MAK,esiM) of getting 
an ellipticity emak using MAK displacement and an ellip- 
ticity esiM) using simulated displacement. It is important 
to check P(£MAK,esiM) if, as we will do in future, we want 
to estimate cosmological parameters from ellipticity distri- 
bution. This function tells us what accuracy we may expect 
on the ellipticity measurements. We represent this probabil- 
ity in the left panel of Figure [S] We see a strong correlation, 
already seen in Fig.[4l indicating a high accurate reconstruc- 
tion. We also see that the error seems low. We represent at 
the left side of the thick red solid line of the right panel the 
conditional probability P(eMAK|£siM) that we computed us- 



ing: 



D/ I \ -P(eMAK, £SIm) 

P(£mak|£sim) = -i — ^ — '— (35) 

_,, P(£,£SIM)d£ 



wherever it was possible to evaluate the denominator. This 
conditional probability is mostly Gaussian, so we tried to 
estimate the standard error by computing the mean variance 
of the error on the ellipticity for the distribution between the 
two dotted red line £sim G [es,min;£s,-m,ax\ with es,-m,in = 
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Figure 4. MAK reconstructed ellipticity distribution us. ellipticity distribution from simulated displacement field - This figure gives tlie 
ellipticity distribution computed using eitlier tlie MAK reconstructed displacement field (solid black line) or tlie simulated displacement 
field (solid thick gray line). We represent the dispersion expected if the error on estimating the distribution come from the number of 
voids in a given bin. We assumed a Poisson distribution for the estimation of the error bar. The displacement fields were both smoothed 
with a Gaussian kernel of radius Rf = 2.5/i~^Mpc in the left panel and to E.f = 4/i~^Mpc in the right panel. 
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Figure 5. Ellipticity derived from the simulated displacement field vs. the MAK displacement field - We represent here a scatter plot 
between the cUipticities of the voids that were both detected in the MAK reconstructed displacement field and the simulated displacement 
field, both smoothed at 4/i~^Mpc. In the left panel, we show the raw joint probability distribution of the two ellipticities. The gray-scale 
is linear in the density of probability. In the right panel, we represent the conditional distribution of Emak given ssimi on the left of the 
thick vertical line. On the right of this same line, we represent this same distribution if one uses the estimated standard deviation o-g of 
the error. It is estimated using the distribution between the two vertical dotted lines. The estimated standard deviation is = 0.018. 
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eoivA 

Figure 6. Tidal ellipticities vs Volume ellipticity - This plot 
gives a comparison of the elhpticity of the void as determined 
either using the shear of the displacement field [eDlVAi Eq. I I13II ] 
or using the overall shape of the void [evoh Eq. IjlSp]. 



0.15 and es,max ~ 0.40. With this notation, we computed 
1 

O-e = 




-e)2P(eMAK|e). (36) 



Within the interval limited by the two dotted red line, we 
estimate ~ 0.018. At the right of the vertical red sohd 
line, we show the function 

P{e\esiMv,(^e) = — e ^'^ (37) 

with as estimated above. We note that this model of 
the conditional probability (right of the vertical solid line) 
looks much like the actual ellipticity discrepancy (left of the 
vertical solid line). 



5.3.3 Volume ellipticity Svoi vs. Tidal ellipticity ediva 

In Fig. |6l we represented a comparison between the elliptic- 
ity of the volume, EvoI, and the local tidal ellipticity, ediva. 
The volume ellipticity is computed using the Eq. psp . with 
the inertial mass tensor M as computed using the smoothed 
displacement field. Visually, the two variables seem loosely 
correlated. We observe they do follow each other but this 
correlation just get worse when the ellipticity increases. This 
is expected as the volume ellipticity is a non-local quan- 
tity and thus is sensitive to all local ellipticities in the void 
volume. This is what makes EvoI more difficult to use as a 
precise cosmology probe. We show in Appendix[C]that these 
two quantities are indeed related but only to first order. This 
explains that the scatter seems smaller for small ellipticities 
than for high ellipticities. The volume ellipticity, which has 



been used so far, seems thus to he a poor proxy of the tidal el- 
lipticity, which we manage to recover with extreme precision 
using our Lagrangian orbit reconstruction technique. For the 
rest of this paper, we will only use the tidal ellipticity. 

5.3.4 Simulation vs Analytic 

In this section, we compare the results obtained on the simu- 
lated displacement field and the prediction given by the an- 
alytical model of Section |4l We represented in the left panel 
of Fig. [5] the ellipticity distribution of all observable voids 
as defined in Section [3. II The black points give the elliptic- 
ity distribution for voids in the reconstructed displacement 
field. We estimated error bars assuming a Poisson distri- 
bution of the number of voids in each bin. The red line is 
obtained using the method of Appendix iBl The dashed blu e 
line is obtained through the formula of lPark fc Led |2003), 
where no constraints are put on the curvature of the density 
field where the ellipticity of the void is measured. 

The success of the comparison between the black points 
and the solid red curve is striking. It shows the importance 
of our constraint that we only look in the negatively curved 
part of S*. We note that this should be a robust feature 
for voids found with any void finder. Any probe of the el- 
lipticity in voids looks in regions of the density field that 
is strongly underdense, and thus should come mainly from 
initially underdense regions. In these primordial underdense 
regions, the curvature of the density field is likely to be pos- 
itive, which corresponds to a negatively curved in our 
case. Our measured ellipticity distribution in the simulation 
IS very clean because we rely on features of the displacement 
field which can be understood in terms of linear theory even 
at redshift 2 = 0. 

In the right panel of Fig. (2] we show this same ellipticity 
distribution but only for "true voids". The comparison be- 
tween simulation and analytic is also here successful. As we 
noted in Section [4.31 there is no real difference between the 
two models in this case. However, there is no way a purely 
geometrical analysis could yield this curve from the obser- 
vation of galaxy catalogues, as this requires the knowledge 
of the sign of eigenvalues of Ti^m (Eq. [6]). We note a small 
shift of ~l-3% between the model and the reconstruction. 
We find, using numerical experiments with Gaussian ran- 
dom fields, that a fraction of this shift may be explained 
by the finite bin size and the very strong steepness of the 
distribution represented in this panel. 

5.4 Evolution with redshift 

In this section, we focus on the evolution with redshift of 
the ellipticity of voids. This evolution h as been sh own to be 
analytically a sensitive probe of w by iLee fc Park (200^. 
We took snapshots of the simulation at different redshifts 
and computed the ellipticity distribution in each of these 
snapshots. We note that we must at least have two main dif- 
ferences compared to what would happen with observations. 
First, we may have a systematic effect in the evolution of the 
mean ellipticity as we are studying only a single realisation 
of initial conditions. Second, the number of available voids 
should be a non-trivial function of redshift. Second, because 
of both volume and selection effects, the error bars should 
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Figure 7. Evolution of the mean ellipticity with redshift - We represent the evolution of the mean ellipticity with redshift for the two 
ACDM simulations (left panel) and the wCDM simulation (right panel) . The mean ellipticities deduced from ASIM are represented using 
square symbols, and the ones from ASIM2 using triangular symbols. The solid curve is obtained using the statistical model of AppendixlB] 
and changing erg according to the growth factor as specified in Eq. I I38I I. The bottom panels give the relative difference, in percentage, 
between the simulation and the analytical model. In the bottom left panel, the points at ~2% correspond to the square symbols 



be large at both small and large redshift, while attaining a 
minimum at some intermediate redshifts The exact numbers 
for this second problem depends on the specific considered 
galaxy survey, in particular the apparent magnitude limit 
and the incompleteness function. 

To compute the predicted ellipticity distribution at any 
given redshift z, we simply scale the 0-8(2:) parameter using 
the growth factor D{z): 

as(z) = as{z = 0)x^^^^. (38) 

For clarity we only represent on Fig. [7] the mean ellip- 
ticity e, defined as, 

e{z) = / eP{e\z)de (39) 
Jo 

with P{e\z) the probability distribution of the ellipticity e 
at redshift z. The red solid line gives the prediction given 
by the analytic model of Section [4] The black points are 
obtained from the simulated field. The error bar on the mean 
ellipticity is estimated using 

(40) 

V 'voids 

with cTg = 0.02 as estimated in Section [5. 3. 2 1 for a smoothing 
scale of 4/1^ ^Mpc. For A^voids — 10"', this gives a typical error 
of CTe = 2 X 10~^ on the mean. This error bar corresponds 
to the uncertainty of the ellipticity derived from the MAK 
reconstructed displacement field with respect to the one ob- 



tained from the simulated displacement field. This gives an 
interval on which the mean ellipticity can be trusted. 

In the left panel of Fig. [T] we see that the agreement 
between the analytical model and the one from the simu- 
lated displacement field (square symbols, "Simulation 1") 
is very good for w — —1. Looking at the relative error be- 
tween the model (lower-left panel) and the simulation yields 
a systematic ~2% deviation relative to the expectation. The 
main reason is the inexactness of the initial conditions to 
the finite number of modes available in the simulation box. 
Even though the power spectrum is normalised to as = 0.77, 
the realized ag of the displacement field is 0.783. This pro- 
duces intrinsically a shift of 1.7% towards positive errors. It 
is exactly what we observe. As we see, this bias is relatively 
modest. However it should be observable when we look the 
small amplitude of the expected reconstruction errors given 
by the small error bars. To check this effect, we rerun an- 
other simulation with the exact same cosmology but with 
another seed. This time, we measured as = 0.7688 in the 
initial conditions, this corresponds to a small statistical fluc- 
tuation of —0.15%. We plotted the corresponding evolution 
of the mean ellipticity in the left panel of Fig. [7] (triangular 
symbols, "Simulation 2"). 

This will not prevent applying this method to obser- 
vations for two reasons. First, we will marginalise over the 
bias and so the systematic shift will disappear. Second, each 
considered slice should be a nearly independent random re- 
alisation of a Gaussian random field normalised to the same 
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Us. Thus the points should be scattered according to our 
dashed horizontal line "0%" and not systematically pushed 
up or down. 



5.5 w = —0.5 vs TO = —1.0 

In all the previous sections, we studied the case of a stan- 
dard ACDM cosmology with w — —1. However, we first 
started to look at voids to check if they may be good tracers 
of the properties of the dark energy, and in particular of the 
growth factor. We now focus on the results obtained from 
wSIM, a wCDM simulation with w = —0.5 as we specified 
in Section [5.11 The results are presented in the right panel 
of Fig. [7] and in Fig. [S] We computed that our particular 
realisation of the initial conditions had ag = 0.773, which is 
0.33% above 0.77. We again note that the discrepancy in the 
lower-right panel in Fig. [S] has the correct systematic shift 
at low redshift. Taking into account this shift, as previously, 
the analytical model seems to follow the results given by 
simulation at the 0.1% level, taking into account the statis- 
tical uncertainty. The points obtained from the simulation 
are in excellent agreement with the simulation. 

Current redshift galaxy catalogues map the Universe at 
intermediate redshifts < z < 1. We check if our method 
is sufficiently sensitive to distinguish a. w = —0.5 from a 
w = —1 cosmology in Figure [S] In this figure, we used the 
erg as measured in the simulations to compute the analytical 
predictions (red and blue solid curves). We note that even 
at z ~ 0.2, the behaviour of the two curves is already sig- 
nificantly different and above statistical uncertainties. If we 
consider the whole interval between 2 = and z = 1, the 
difference is very significant compared to the uncertainties, 
with an ellipticity that changes by ~35%. 

In all the above we considered catalogues with an im- 
portant observable number of voids (typically ~10,000). We 
do not expect to have such a high number available in cata- 
logues. Now we try to make an estimate of the error bars on 
the mean ellipticity that we may expect. The SDSS covers 
one fifth of the sky. We now limit the survey at z = 0.1 
(~300/i"^Mpc), and we take a Lagrangian smoothing scale 
of 5/i~^Mpc. This smoothing scale is motivated by the av- 
erage density of galaxies i n the SDSS, which i n band r is 
about 1 - 5 X 10"^/i^Mpc^ (|Blanton et al.ll2003l '). This gives 
a mean separation of ~ 2 — 4/i~^Mpc. The equation (jSip 
predicts that we should observe ~1,000 of our voids in this 
volume when smoothing the density field in Lagrangian co- 
ordinates with a Gaussian of radius 5/i~^Mpc, taking into 
account the survey coverage. If we go to z = 0.2, this number 
should increase to ~9,000. This means that the error bars 
should only be moderately larger than the one that we con- 
sidered in this work. Typically we expect about three times 
larger. Even, with this amount of uncertainty, the compari- 
son to the analytic model should be able to highly constrain 
the equation of state of dark energy at z < 0.2. The Fisher 
matrix analysis is done in a companion paper. 



6 COMPARISON OF DIVA TO EARLIER 
VOID FINDERS 

In this section, we discuss how our technique is related to 
the other existing void finders. We try to make a qualitative 
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Figure 8. Difference between w=-l and w=-0.5 - We plotted 
here the evolution of the mean ellipticity with redshift. We used 
the simulation ASIM2 (square) and wSIM (triangle). The red solid 
line gives the prediction for erg = 0.7688, w = —1. The blue solid 
line gives the prediction erg = 0.773 and w = —0.5. These two 
values of erg have been computed using the initial conditions of 
the two simulations. 



assessment of its strengths and weaknesses compared to the 
three classes of void finders define in Section [T] 

The void finders of the first class try to find emp- 
tier regions in a distribution of points, which in ac- 
tual catalogue s correspond to ga laxies. The void finder 
developed b y lEl-Ad et al.l (Il997f), and one of its later 
versions by iHoyle fc Vogelev j2C)02h. i s popularly used 
in observations (iHovle fc Vogeleyl 12004 iHovle et al.l [ioosi : 
iTikhonov fc Karachentsevll2006l : [Foster fc Nelsonll2009l ). In 



these void finders, the first step consists in classifying galax- 
ies in two types. Galaxies may lie in a strongly overdense 
region, in this case we consider it as a "wall galaxy". The 
other possibility is that they lie in a mildly underdense re- 
gion, and they are then called "field galaxies". The exact 
separation between "wall galaxies" and "field galaxies" de- 
pends on an ad-hoc parameter. This parameter specifies how 
the local density of galaxy control the type (field or wall) of 
the galaxy. The voids are then grown from regions empty 
of wall galaxies. This classification thus gives a non-trivial 
dependence of the void sizes and shapes on the galaxy bias 
and catalogue selection function. Additionally, this defini- 
tion depends on an ad-hoc empirical factor. These issues 
make the quantitative study of the geometry of voids dif- 
ficult, while they find voids that correspond to the visual 
impression given by large scale structures in redshift cata- 
logues of galaxies. 

Void finders belonging to the second class look for par- 
ticularities in the continuous three dimensional distribution 
of the dark matter traced by galaxies. Of course, from obser- 
vational data, one has then to first project and then smooth 
the distribution of galaxies to obtain this distribution. Dif- 
ferent techniques are used; 
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- One technique is to adaptively smooth the galaxy 
distribution either using an SPH technique (see e.g. 
IColombi et"aDl2007l) or a Delaunay Tesselation Field Esti- 
mator (jSchaap &: van de Wevgaer'?l200(]| ). Voids must then 



be identified from the smooth distribution derived using ei- 
ther of these techniques. O ne option is to use a sche me simi- 
lar to lEl-Ad et all (|l997l ) or lHovle fc VogelevI (|2002l ) to iden- 
tify shapes of underdensity in the vicinity of a minima of the 
density field (|Colberg et all 120051 '). As with void finders of 
the first class, the galaxy bias does not affect the positions 
of the voids but their overall pr operties. A second option is 
to use a Watershed Transform (|Platen et al.|[2007l ) to iden- 
tify the cavities of the voids. In this case, the galaxy bias 
does not affect the structure of the cavities. However, devis- 
ing an efficient way of relating the geometrical properties of 
these cavities to the cosmology, which corresponds to study- 
ing Morse theory, could well be non-trivial (, Jgst 20Q^) . Some 
work to study the skeleton (also called the "cosmic spine") 
of Large Scale structures in this t h eory has been r e cently 
don e bvlAragon-Calvo et all (|2008l '). iPogosvan et all Hooi) 
and lSousbie et all (|2009l v 

- A second technique consists in using the local density 
estimated using the Vor onoi diagram of a Delaunay tessel- 
lation to locate minima (|Nevrinckll200l ). The particles are 
first grouped in zones. Each particle is assigned to a zone on 
to where this particle is attracted if it followed the density 
gradient as in the watershed technique. Each zone is defined 
to be a void. But it is also possible to define a hierarchy 
of voids by checking, for two neighbouring voids, which of 
the two has the lowest density at the minima. The zones are 
then grouped and a probability of being a void is assigned 
depending on the contrast between the density at the ridge 
of the void and its depth. 

This void finder has the advantage of defining voids in terms 
of topology of the density field, which is easier to handle from 
a theoretical point of view and may better define a void 
in terms of dark matter. Still, we are faced with the task 
of relating the shape of the voids that are found by these 
algorithms, which is non-local by nature, to cosmology. As 
mentioned previously for void finders of the first class, this 
seems to be non-trivial. 

Void finders of the third class use the inferred dark mat- 
ter distribution but they do a dynamical analysis to infer the 
location of these voids. DIVA belongs to these class of void 
finders. There are two advantages of looking at dynamics 
for voids, (i) It gives a much more physical and intuitive 
definition of these structures: voids corresponds to places in 
the universe from which the matter is really escaping and 
not gravitationally unstable at present time, (ii) Using this 
criterion, one is bound to use either the velocity field or the 
displacement field. These two quantities are h ighly linear. It 
has be en directly shown for vel ocity fields by Ciecielg et al.l 
(|2003|) and indirectl y shown bv lMohavaee et al.l ( 20061 ) and 
iLavaux et al.l (|2008l ) for the displacement field. This linear- 
ity helps us at constructi n g an analy tical statistical model 
of the voids. iHahn et al.l ( 20071 ) and iForero- Romero et al.l 
(|2009l ) attempted to classify structures according to a cri- 
terion on the gravitational field. However we may highlight 
two very important differences compared to the approach 
we are following here: 

- we are using a purely Lagrangian method and it takes 



into account the true evolution of the void and not how 
virtual tracers in the void should move now, 

- we put an exact natural thresh old on eigenvalues to 
classify voids. This is in contrast with lForero-Romero et ali 
(|2009h where they need to to put a threshold on eigenval- 
ues depending on an estimated collapse time. In our case, 
everything is already integrated in the definition of the dis- 
placement field. 

Moreover, the Monge-Ampere-Kantorovitch reconstruction 
presents the two advantages of: (i) never diverging in the 
neighbourhood of large density fluctuations, compared to 
a pure gravitational approach, (ii) recovering exactly the 
Zel'dovich approximation in the neighbourhood of centres 
of voids. 

As for the other void finders, DIVA depends on galaxy 
bias. We recall that the bias b is defined with 



5g ~ bSm, 



(41) 



with Sg the density fluctuations of galaxies and 5,n the den- 
sity fluctuations of matter. As MAK is essentially recon- 
structing the Zel'dovich displacement in underdense regions, 
and that the Zel'dovich potential is proportional to density 
fluctuations, the MAK displacement should also be mostly 
linear with the bias. 

We describe in a second paper (Lavaux & Wandelt 2009, 
in prep.) how to relate the volume of the voids that we flnd 
in Lagrangian coordinates to the voids that we observe in 
Eulerian coordinates. 



7 CONCLUSION 

We have described a new technique to identify and char- 
acterise voids in Large Scale structures. Using the MAK 
reconstruction, we have been able to define void centres in 
Lagrangian coordinates by assimilated them to the maxima 
of the divergence S<s of the displacement field, interpreted 
as its source term. The scalar field has the interesting 
property of being nearly equal to the opp osite of the lin- 
early extrapolated primordial density field (|Mohavaee et al.l 
120061 ). This allowed us to consider that the statistics of those 
two fields were equal. Using that, we made predictions on 
the number of voids in Lagrangian coordinates, along with 
their ellipticities defined using the eigenvalues of the curva- 
ture of . 

We tested our model using A'"-body simulations with 
different cosmologies {w = —1 and w = —0.5). We checked, 
using the largest Lagrangian reconstruction run so far, that 
MAK is capable of recovering the ellipticity of individual 
voids to the order of a few percent. We highlighted the im- 
portance of using our model for the statistics o f the eigenval- 
ues of the curvature instead of the formula of iDoroshkevichI 
1 19701 ) for the particular case of voids. We showed that our 
analytical model agrees within ~0. 1-0.3% to results obtained 
with MAK and the displacement field obtained from the sim- 
ulation. We expect our method to be able to provide a very 
promising constraint on the equation of state of the Dark 
Energy of the late universe, especially given the notable ac- 
curacy of the prediction that we obtained Fig. [S] 

We intend to pursue this work to continue character- 
ising analytically the voids found by DIVA, in particular 
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the evolution of the number of voids and their size distribu- 
tion. We will make further robustness tests using mock cat- 
alogues, including especially redshift distortion effects. We 
also would like to apply our met hod to the Luminous R ed 
Galaxy sample of the SDSS DR7 (|Abazaiian et al.ll2009l ). 



ACKNOWLEDGEMENTS 

This research was supported in part by the National Sci- 
ence Foundation through TeraGrid resources provided by 
the NCSA. TeraGrid systems are hosted by Indiana Uni- 
versity, LONI, NCAR, NCSA, NICS, ORNL, PSC, Purdue 
University, SDSC, TACC and UC/ANL. The authors ac- 
knowledge financial support from NSF grant AST 07-08849. 
We acknowledge the hospitality of the California Institute of 
Technology, where the authors completed most of this work. 
We thank the referee, Rien van de Weygaert, for his useful 
comments and suggestions. 



REFERENCES 

Abazajian K. N., Adelman-McCarthy J. K., Agiieros M. A., 

et al. 2009, ApJS, 182, 543 
Aragon-Calvo M. A., Platen E., van de Weygaert R., Szalay 

A. S., 2008, ArXiv e-prints 
Bardeen J. M., Bond J. R., Kaiser N., Szalay A. S., 1986, 

ApJ, 304, 15 

Benoist C, Maurogordato S., da Costa L. N., Cappi A., 

Schaeffer R., 1996, ApJ, 472, 452 
Benson A. J., Hoyle F., Torres F., Vogeley M. S., 2003, 

MNRAS, 340, 160 
Bertschinger E., 1983, ApJ, 268, 17 
Bertschinger E., 1985, ApJS, 58, 1 
Bertschinger E., 2001, ApJS, 137, 1 

Bertsekas D. P., 1979, A Distributed Algorithm for the As- 
signment Problem. MIT Press, Cambridge, MA 

Blanton M. R., Hogg D. W., BahcaU N. A., Brinkmann al. 
2003, ApJ, 592, 819 

Blumenthal G. R., da Costa L. N., Goldwirth D. S., Lecar 
M., Piran T., 1992, ApJ, 388, 234 

Bond J. R., Myers S. T., 1996, ApJS, 103, 1 

Brenier Y., Frisch U., Henon M., Loeper G., Matarrese S., 
Mohayaee R., Sobolevskii A., 2003, MNRAS, 346, 501 

Catlett C, Hallcock W. E., Andrews P., Aydt R., al. 2008, 
in Lucio Grandinetti ed., HPC and Grids in Action Vol. 16 
of Advances in Parallel Computing, TeraGrid: Analysis of 
Organization, System Architecture, and Middleware En- 
abling New Types of Applications, pp 225-249 

Ciecielg P., Chodorowski M. J., Kiraga M., Strauss M. A., 
Kudlicki A., Bouchet F. R., 2003, MNRAS, 339, 641 

Colberg J. M., Pearce F., Foster C, Platen E., al. 2008, 
MNRAS, 387, 933 

Colberg J. M., Sheth R. K., Diaferio A., Gao L., Yoshida 
N., 2005, MNRAS, 360, 216 

CoUess M., Dalton G., Maddox S., Sutherland W., Norberg 
P., al. 2001, MNRAS, 328, 1039 

Colombi S., Chodorowski M. J., Teyssier R., 2007, MN- 
RAS, 375, 348 

Croft R. A. C, Gaztanaga E., 1997, Mon. Not. Roy. Astron. 
Soc, 285, 793 



Davis M., Efstathiou G., Frenk C. S., White S. D. M., 1985, 
ApJ, 292, 371 

de Lapparent V., Geller M. J., Huchra J. P., 1986, ApJL, 
302, LI 

Doroshkevich A. G., 1970, Astrophysics, 6, 320 
Efstathiou G., Frenk C. S., White S. D. M., Davis M., 1988, 

MNRAS, 235, 715 
El-Ad H., Phan T., Dacosta L. N., 1997, MNRAS, 287, 790 
Erdogdu P., Huchra J. P., Lahav O., et al. 2005, MNRAS 
Forero- Romero J. E., Hoffman Y., Gottlober S., Klypin A., 

Yepes G., 2009, MNRAS, 396, 1815 
Foster C, Nelson L. A., 2009, ApJ, 699, 1252 
Goldwirth D. S., da Costa L. N., van de Weygaert R., 1995, 
MNRAS, 275, 1185 
Gottlober S., Lokas E. L., Klypin A., Hoffman Y., 2003, 
MNRAS, 344, 715 
Gregory S. A., Thompson L. A., 1978, ApJ, 222, 784 
Hahn O., Porciani C, CaroUo C. M., Dekel A., 2007, MN- 
RAS, 375, 489 

Hausman M. A., Olson D. W., Roth B. D., 1983, ApJ, 270, 
351 

Hoffman G. L., Salpeter E. E., Wasserman I., 1983, ApJ, 
268, 527 

Hoffman Y., Shaham J., 1982, ApJL, 262, L23 

Hoyle F., Rojas R. R., Vogeley M. S., Brinkmann J., 2005, 

ApJ, 620, 618 
Hoyle F., Vogeley M. S., 2002, ApJ, 566, 641 
Hoyle F., Vogeley M. S., 2004, ApJ, 607, 751 
Huchra J. P., 2000, in Courteau S., Willick J., eds. Cosmic 

Flows Workshop Vol. 201 of Astronomical Society of the 

Pacific Conference Series, The 2MASS Redshift Survey. 

pp 96-+ 
Icke v., 1984, MNRAS, 206, IP 

Joeveer M., Einasto J., Tago E., 1978, MNRAS, 185, 357 
Jones D. H., Read M. A., Saunders W., et al. 2009, MN- 
RAS, 399, 683 

Jost J., 2008, Riemannian Geometry and Geometric Anal- 
ysis. Universitext, Springer Berlin Heidelberg 

Kauffmann G., Fairall A. P., 1991, MNRAS, 248, 313 

Kirshner R. P., Oemler Jr. A., Schechter P. L., Shectman 
S. A., 1981, ApJL, 248, L57 

Lavaux G., Mohayaee R., Colombi S., TuUy R. B., 
Bernardeau F., Silk J., 2008, MNRAS, 383, 1292 

Lavaux G., TuUy R. B., Mohayaee R., Colombi S., 2008, 
ArXiv e-prints 

Lee J., Park D., 2006, ApJ, 652, 1 

Lee J., Park D., 2009, ApJL, 696, LIO 

Lewis A., ChaUinor A., Lasenby A., 2000, ApJ, 538, 473 

Martel H., Wasserman I., 1990, ApJ, 348, 1 

Mohayaee R., Mathis H., Colombi S., Silk J., 2006, MN- 
RAS, 365, 939 

Neyrinck M. C, 2008, MNRAS, 386, 2101 

Norberg P., Baugh C. M., Hawkins E., al. 2001, MNRAS, 
328, 64 

Novikov D., Colombi S., Dore O., 2006, MNRAS, 366, 1201 
Park D., Lee J., 2007, Physical Review Letters, 98, 
081301-h 

Park D., Lee J., 2009a, MNRAS, pp 1323— h 

Park D., Lee J., 2009b, MNRAS, 397, 2163 

Patiri S. G., Betancort-Rijo J. E., Prada F., Klypin A., 

Gottlober S., 2006, MNRAS, 369, 335 
Peebles P. J. E., 1989, ApJL, 344, L53 



Precision cosmology with voids: definition, methods, dynamics 17 



Percival W. J., Nichol R. C, Eisenstein D. J., et al. 2007, 
ApJ, 657, 645 

Piran T., Lecar M., Goldwirth D. S., da Costa L. N., Blu- 
menthal G. R., 1993, MNRAS, 265, 681 

Platen E., van de Weygaert R., Jones B. J. T., 2007, MN- 
RAS, 380, 551 

Platen E., van de Weygaert R., Jones B. J. T., 2008, MN- 
RAS, 387, 128 

Plionis M., Basilakos S., 2002, MNRAS, 330, 399 

Pogosyan D., Pichon C, Gay C., Prunet S., Cardoso J. F., 
Sousbie T., Colombi S., 2009, MNRAS, pp 673— h 

Regos E., Geller M. J., 1991, ApJ, 377, 14 

Ryden B. S., 1995, ApJ, 452, 25 

Ryden B. S., Melott A. L., 1996, ApJ, 470, 160 

Sahni V., Sathyaprakah B. S., Shandarin S. P., 1994, ApJ, 
431, 20 

Sahni V., Shandarin S., 1996, MNRAS, 282, 641 

Schaap W. E., van de Weygaert R., 2000, A&A, 363, L29 

Shandarin S., Feldman H. A., Heitmann K., Habib S., 2006, 

MNRAS, 367, 1629 
Sheth R. K., van de Weygaert R., 2004, MNRAS, 350, 517 
Sousbie T., Colombi S., Pichon C, 2009, MNRAS, 393, 457 
Tegmark M., Blanton M. R., Strauss M. A., al. 2004, ApJ, 

606, 702 

Tegmark M., Eisenstein D. J., Strauss M. A., Weinberg 

D. H., Blanton M. R., al. 2006, PhRvD, 74, 123507 
Teyssier R., 2002, A&A, 385, 337 

Tikhonov A. V., Karachentsev I. D., 2006, ApJ, 653, 969 
TuUy R. B., Fisher J. R., 1978, in Longair M. S., Einasto J., 
eds. Large Scale Structures in the Universe Vol. 79 of lAU 
Symposium, Nearby small groups of galaxies, pp 31-45 
van de Weygaert R., van Kampen E., 1993, MNRAS, 263, 
481 

York D. G., Adelman J., Anderson Jr. J. E., Anderson 

S. F., Annis J., al. 2000, AJ, 120, 1579 
Zel'dovich Y. B., 1970, A&A, 5, 84 



APPENDIX A: EIGENVALUES PROBABILITY 
DISTRIBUTION 

In this appendix, we derive the unconditional distribution of 
the eigenvalues of the Jacobian matrix of the displacement 
field in the Zel'dovich regime. Starting from the potential 
$(q) of the displacement field at the comoving Lagrangian 
coordinate q, we define the Jacobian matrix of the displace- 
ment field 



the matrix V of the variances 



Ti, 



dq,dqj 



(Al) 



This matrix is real and symmetric. We assume the com- 

J ionents to be norma lly distributed and to be as given by 
Bardeen et al.lll986l . Appendix A): 
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thus the covariance matrix C is 



C = V~ 
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(A4) 



The unconditional probability dP{{TA}) of observing these 
components at a point in the universe is given by 

dPi{TA}W) = J^^l exp f-i fl Ca,bTaTb] fl dTA 
16V5a67r3 \^ 2 J 

(A5) 

We want now to make a change of variables and get the prob- 
ability distribution of the eigenvalues of Tij. The appendix 
B ofyBardccn ct al. (198^) shows that the infinitesimal vol- 
ume is transformed according to: 

6 

H dTA = l(Ai - A2)(Ai - A3)(A2 - \b)\ dAidA2dA3dns3 , (A6) 

A^l 

with di}g3 the infinitesimal rotation of the hypersphere of 
dimension 3. The quadratic form Q in the argument of the 
exponential of Eq. (IA5P may be expanded 



1 5 

= Ca,bTaTb = — {^{T, 



-To" 



T|) 



-2(riT2 + TiT-i + T2T3) + r| + t| + T^) 
Linear algebra tells us that the scalar quantities 

Ki = T1+T2+TS and K2 



(A7) 



TiT2+T2T3+TiT3-T'i-Ti-Ti 

(A8) 

are invariant by change of vector basis. The quadratic form 
may then be rewritten: 



Q 



(A9) 



We may now express Ki and K2 in terms of the eigenvalues: 

Ki = A1+A2 + A3, (AlO) 
K2 = A1A3 + A1A2 + A2A3. (All) 



After integrating the expression (|A5|l over 5* , we obtain the 
unconditional probability that the Jacobian matrix has the 
three ordered eigenvalues Ai, A2 and A3: 



(A2) 



with a the standard deviation of the density field. As it 
is a 3 X 3 real symmetric matrix, there are only 6 in- 
dependent components. We label these components with 
A = 1 . . . 6, where each number refer to the couples 
(1, 1), (2, 2), (3, 3), (1, 2), (1, 3) and (2, 3). We may now write 



P(Ai,A2,A3ia) 



3375 



exp 



3 (2Kt - 5K2) 



2(72 



|(Ai-A2)(Ai-A3)(A2-A3)|. 

(A12) 

This equation corre sponds to the distribution derived by 
iDoroshkeviciil (| 19701 ). 
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APPENDIX B: PROBABILITY DISTRIBUTION 
OF TIDAL FIELD IN VOIDS 

We show in this Appendix how to compute the form of the 
probability distribution of the gravitational tidal field in the 
case of voids. The technique involved in this derivation looks 
much like the one used in Appendix |X] but with an extra 
complication due to correlations with the density field. By 
definition, void centres are chosen to be maxima of 



(Bl) 



with $ defined as in § [2l the scalar potential of the displace- 
ment field. We will assume that, for voids, we are fully in 
the Zel'dovich regime and thus we can equate $ as given 
by MAK to the potential of the Zel'dovich displacement at 
high redshift. In this case, 5* corresponds the primordial 
density fluctuations scaled linearly to 2 = 0. 

Being a maxima of means the displacement field 
satisfies two properties: the gradient of is null and the 
Hessian matrix 

as* 



H 



dqidqj 



(B2) 



of S<s is negative-definite. Thus our aim is to compute the 
probability of the Jacobian matrix of the displacement field 
given that Hi^j is symmetric negative-definite. We write this 
probability P{Ti,j\Hi,^ < 0). 

We assume that the gravitational potential is a Gaus- 
sian random field determined by the power spectrum of mat- 
ter density fluctuation P{k), with k a wave number. We 
compute the correlation between those two flelds: 



with 



^2 

— {SijSl^rn + 5i,l5i,m + 5i,m5j,l) , (B3) 



— {5i^j5l_rn + Si^lSj.m + Si.mSj.l) , (B4) 



(SijSl^m + Si,lSj^m + Si^mSj,l) , (B5) 



2 



S2 = <J^ and ajj = 56 and Tth 



-S4, 



a the variance of the density fleld and 
1 



IOtt^ 



k''P{k)dk. 



(B6) 



(B7) 



To reduce the complexity of the correlation matrix, we now 
use the reduced random variables defined as follow 

T — -^''J 



Hi 



(B8) 
(B9) 

(BIO) 



and the reduced correlation 

r — . 

arUH 

As T and H are 3x3 real symmetric matrices, there 
are only 6 independent components for each of these matri- 
ces. As in the Appendix]^ we label these components with 



A = 1 . . . 6, where each number refer to the couples 
(1,1), (2, 2), (3, 3), (1,2), (1,3) and (2,3). The matrix V of 
the variance of the 12 reduced components may be formally 
written, using a block-matrix representation: 



V = 



and 



A 



A rA 
rA A 



(BU) 
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\ 1 / 



(B12) 



The inverse may be computed straightforwardly 
1 



C = = 



1 



A-' 
-rA-'' 



-r-A-^ 
A-^ 



(B13) 



Now, we may express the joint probability P{T, H) of 
observing a tidal field T for the gravitational potential and a 
curvature H for the density field using the covariance matrix 
C = M-^. 



P{T,H\r)^ 



v/|dctCl 
(27r)6 



-exp 



(B14) 



with Y — {T,H). Now the probability of observing some 
matrix T given that H must be positive could be computed 
formally: 



P(r|r, H <Q)= [ P{T, H\r)AH 

J H<0 



(B15) 



However, it is quite involved to find an analytic expression 
of this integral in function of the eigenvalues of T. We pro- 
pose to sample this distribution instead of computing this 
integral. 

One can prove that the conditional probability P{T\H) 
may be written 



P{T\H,r) 



PiT,H\r) 
S^P{T,H\r) 



exp 



\T-rH)A{T-rH)] (B16) 



(27r)3 \ 2 

Thus, Equation (|B15|I may be re-expressed as: 



P(T\r,H < 0) = 



P{T\r,H)P{H)dH 



(BIT) 



with P{H) the probability of getting a random symmetric 
matrix H, with the covariance matrix A. The method is now 
the following: 

- we generate a random symmetric matrix H , if it is neg- 
ative, we accept it, in the other case we try again; 

- we generate a matrix T, following the probability given 
by Equation (|B16P : 

- we compute and store the eigenvalue of this matrix T, 
after multiplication by ar- 

That way the joint probability distribution 
P{Xi, \2, ^a\(yT ,r, H < 0) is correctly sampled, even 
though we do not have any explicit expression of it. 



Precision cosmology with voids: definition, methods, dynamics 



APPENDIX C: LOCAL TIDAL ELLIPTICITY 
VS. GLOBAL VOLUME ELLIPTICITY 

In this appendix, we try to relate the two elhpticities EvoI 
and ediva defined in Eq. (|13|) and (|15|) . To do that, we will 
make use of Zel'dovich approximation in voids, which has 
been shown to be a relatively precise modeling of the void 
evolution. The inertial mass tensor writes as: 

M = aT-K (CI) 

with I the 3x3 identity matrix, 

d^q||x(q)-x||^ (C2) 

IV 

i'^q(a::,(q) - xO(Kj(q) - Sj), (C3) 

IV 

with q the Lagrangian coordinates, V the Lagrangian do- 
main of the considered void, x the centre of mass of the 
V in Eulerian coordinates. With this parametrization, the 
volume ellipticity Svoi simplifies as 

/ T \ 1/4 

e..i = l-r-j-\ (C4) 

with Ji and J3 the smallest and largest eigenvalues of K. 
We may write exactly: 

x(q) = q+*(q). (C5) 

We now expand ^ to first order around the position of the 
centre of mass in Lagrangian coordinates q 

*,(q) = *,(q) + ^fe-qj). (C6) 
aqj 

Using Zel'dovich approximation we identify d'^i/dq^ and 
Ti,j given in Eq. ((8|. We now reexpress Ki^j 

Kij — Lij + Ti^kLk.j + Tj.kLk,i (C7) 

d'q(g.-gO (C8) 

As voids, in Lagrangian coordinates, should be mostly 
isotropic the Lagrangian inertial tensor Lij must be diago- 
nal: 

L^,j = ^aiSij, (C9) 
with 

al^ f d^qllq-qll^ (CIO) 
Jv 

This assumption is verified in average by linear theory but 
may be broken for some specific voids. In the case where 
voids are effectively isotropic, the inertial mass tensor Kij 
is extremely simplified: 

The eigenvalues of K are thus 

J, ocl-f 2A, ~ (1 + AO^ (C12) 

The volume ellipticity may thus be related to the tidal el- 
lipticity as 



